proc buildGroupList {targetAtomPdb selections args} {
    global groups masses
    set targets {}
    set masses {}
    set inStream [open $targetAtomPdb r]
    set index 0
    set targetPdbTop(segnames) {}
    set targetPdbTop(residues) {}
    set idx 0
    set groupby residue
    if {[set idx [lsearch $args by ]] >= 0} {
        set groupby [lindex $args [expr $idx +1] end ]
    }
    switch $groupby {
        segment {
            set group "segment=\$segname"
        }
        residue {
            set group "segment=\$segname,resid=\$resid"
        }
    }
    set index 0
    foreach line [split [read $inStream] \n] {
        set type [string trim [string range $line 0 5]]
        set name [string trim [string range $line 12 15]]
        set resid [string trim [string range $line 22 25]]
        set beta [string trim [string range $line 60 65]]
        set occupancy [string trim [string range $line 54 59]]
        set segname [string trim [string range $line 72 75]]
        if { [string compare $type ATOM ] != 0}  { continue}
        set mygroup [subst $group]
        if ![info exists MyGroups($mygroup)] {
            set MyGroups($mygroup) {}
        }
        foreach selection $selections {
            if [string match $selection $name ] {
                lappend MyGroups($mygroup) $index 
            }
        }
        incr index
    }
    foreach name [lsort [array names MyGroups]] {
        if {[llength $MyGroups($name)] > 0 } {
            lappend groups $MyGroups($name)
        }
    }
    return 
}

proc setupConstraint { forceconstant r0 k0} {
    loadmasses masses
    loadcoords coords
    set zave 0.0
    set zstd 0.0 
    foreach group $groups {
        set comsum {0.0 0.0 0.0 }
        set totalmass 0.0
        foreach atom $group {
        set tmp [vecscale $masses($atom) $coords($atom)]
        set comsum [vecadd $comsum $tmp]
        set totalmass [expr $totalmass +$masses($atom)]
        }
        set com [vecscale [expr 1.0/$totalmass] $comsum]
        set z [lindex $com 2]
        set zave [expr $zave+ abs($z)]
        set zstd [expr $zstd+ $z*$z]
    }
    set zave [expr $zave/[llength $groups]]
    set zstd [expr $zstd/[llength $groups]]
    set zstd [expr sqrt($zstd-$zave*$zave)]
    geoCreate geo ToroidalPore $r0 $k0 [expr $zave-$zstd]
    setGeoRange geo -40.0 40.0 -40.0 40.0 -40.0 40.0
    genGeoConstraint geo Gaussian $forceconstant $zstd $zstd 160
}

